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ABSTRACT 

Energy dissipation is highly intermittent in turbulent plasmas, being localized 
in coherent structures such as current sheets. The statistical analysis of spatial 
dissipative structures is an effective approach to studying turbulence. In this paper, we 
generalize this methodology to investigate four-dimensional spatiotemporal structures, 
i.e., dissipative processes representing sets of interacting coherent structures, which 
correspond to flares in astrophysical systems. We develop methods for identifying and 
characterizing these processes, and then perform a statistical analysis of dissipative 
processes in numerical simulations of driven magnetohydrodynamic turbulence. We 
find that processes are often highly complex, long-lived, and weakly asymmetric in time. 
They exhibit robust power-law probability distributions and scaling relations, including 
a distribution of dissipated energy with power-law index near —1.75, indicating that 
intense dissipative events dominate the overall energy dissipation. We compare our 
results with the previously observed statistical properties of solar flares. 
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1. Introduction 


Intermittency, the inherent inhomogeneity of turbulence, causes energy dissipation to be highly 
localized in space and in time. Spatial intermittency is manifest by coherent structures, while 
temporal intermittency is characterized by irregular, bursty events. Both of these phenomena are 
related, play key roles in the turbulent dynamics, and have important observational consequences. 
As with any dynamical physics problem, a complete understanding of intermittency can only be 
claimed when the solution is described in both the spatial and the temporal dimensions. 

Intermittency is a challenging theoretical problem that impedes our progress toward a complete 
theory of turbulence. Many of the methods used to study turbulence, such as energy spectra and 
correlation functions, are insensitive to intermittency. The methods conventionally employed to 


measure intermittency include structure functions (Muller et al. 2003 Chen et al. 2011; Podesta 


2011), scale-dependent kurtosis ( Wan et ar]|2012 ), topological methods ( Servidio et al.||2009 2010), 
and statistics of discontinuities (Greco et al. 2009a|[b Zhdankin et al. 2012); however, these often 


give incomplete information or are difficult to measure accurately (Jimenez 2007; De Wit 2004 


Sreenivasan fc Antonia||1997 ). Therefore, in order to better guide theoretical models of turbulence, 
it is essential to develop new tools for studying intermittency in a robust and informative manner. 

The statistical analysis of coherent structures is a possible route forward in this direction. It is 
convenient, for both practical and theoretical purposes, to treat coherent structures as discrete ob¬ 
jects due to their localized nature and their central role in the dynamics. The statistical properties 
of these structures, including their intensities and morphologies, give insight into the underlying 
dynamics from which they formed. This information reveals the anisotropy, inhomogeneity, and 
characteristic scales of the dynamics. Furthermore, such an analysis can be extended beyond sim¬ 
ulations and theory, being well suited for a variety of experimental and observational applications. 


including solar flares, the solar photosphere (Cattaneo 1999 Bushby &: Houghton 2005; Stein & 


Nordlund 2006), the interstellar medium (Pan et al. 2009 Falgarone et al. 2015 Boldyrev et al. 


2002; Kritsuk et al. 2011), instabilities in fusion devices (Carbone et al. 2000; Antar et al. 2003 


D’Ippolito et al. 2004), and radiative signatures in optically thin astrophysical plasmas, e.g., in 


black-hole accretion disk coronae (Di Matteo et alT]|1999), hot accretion flows (Eckart et al. 2009), 


jets (Albert et al. 2007), pulsar wind nebulae (e.g. Tavani et al. 2011 Abdo et al. 2011), and hot 
gas in galaxy clusters. 
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The statistical analysis of spatial structures has been a fruitful approach to studying the 
intermittency of turbulence in numerical simulations. The methods for identification of structures 
vary between different studies, but one of the simplest and most common approaches is to identify 
structures as regions in space bounded by an isosurface of the relevant field. For example, this 


was applied in hydrodynamic turbulence to study dissipative vorticity filaments (Jimenez et al. 


1993; Moisy & Jimenez 2004 Leung et al.||2012), which revealed that the radii of intense filaments 


scale with the Kolmogorov microscale while lengths occupy large scales. More recently, similar 
methods were applied to study current sheets and other dissipative structures in simulations of 


plasma turbulence, notably in 2D magnetohydrodynamic (MHD) turbulence (Servidio et al. 2009 


2010 

), 3D MHD turbulence (Uritsky et al.|2010 

Zhdankin et al.|2013 

2014 

), the kinematic dynamo 

(]Wilkin et al. 

2007 

), ambipolar diffusion MHD 

(Momferratos et al. 

2014) 

, boundary-driven MHD 


These previous studies succeeded in describing the morphology and scaling properties of in¬ 
termittent structures, but are incomplete in the sense that they give no information on how the 
structures evolve in time. For this reason, very little can be said about the dynamics, even if 
one assumes statistical stationarity. Important temporal information about intermittent structures 
includes their characteristic timescales, stability, motion, interactions, and impulsiveness. If one 
is interested in understanding these temporal aspects of intermittency, then a broader framework 
must be developed. 


Temporal intermittency was investigated to a limited extent in hydrodynamic turbulence. 
The temporal statistics of vortices in 2D turbulence (Carnevale et al. I991| Whitcher et al. 2008 


Pasquero et al. 2002) and decaying Charney-isotropic geostrophic turbulence (McWilliams et al. 


1999) gave insight into global changes in the population size and morphology of vortices. However, 


the structures were not entirely treated as temporal objects in these studies, which requires an 
algorithm for tracking the structures through time. This can be much more difficult to design and 
implement than the algorithms used to identify and characterize structures in individual hxed-time 
snapshots (Storlie et al. 2004, 2009). To the best of our knowledge, apart from our precursor 


paper (Zhdankin et al. 2015), no systematic studies of spatiotemporal structures were previously 


undertaken for hydrodynamic or MHD turbulence in 3D, and only limited studies were performed 


Aubry et al. 

1991 

Daviaud et al. 

1990 

Jung et al. 

2000 

Golovas & 


Andereck 1997) 


The objective of this paper is to extend the framework previously used for the statistical 


analysis of spatial dissipative structures in numerical simulations of MHD turbulence (Zhdankin 


et al. 2013, 2014) into the temporal realm. We therefore consider the statistical properties of 


4D spatiotemporal objects which represent structures evolving through time. These objects are 
dissipative processes, analogous to flares in astrophysical systems, which, in general, may involve 
many interacting coherent structures. We consider the distributions and scaling relations for process 
characteristics including the total dissipated energy, peak energy dissipation rate, duration, and 
geometric scales; measures of complexity such as the number of interacting structures and types of 
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interactions (mergers and divisions); and the temporal evolution of individual processes. This novel 
methodology is applied to study intermittency in numerical simulations of driven incompressible 
MHD turbulence; however, it can also be applied to hydrodynamic turbulence, kinetic plasma 
turbulence, and other complex dynamical systems exhibiting spatiotemporal self-organization. 


The analysis presented in this paper addresses several fundamental questions regarding the 
intermittency of MHD turbulence. One key question is whether, in the limit of high Reynolds 
number, the overall energy dissipation is dominated by a few intense, long-lasting events residing 
at large scales, or by many weak, short-lived events residing near the dissipation scale. Another 
question is whether there is an inherent relationship between spatial intermittency and temporal 
intermittency, e.g., whether larger structures better retain their coherency in time. A third question 
is whether the dissipative events show any characteristic temporal asymmetry, e.g., impulsive onset 


followed by a slow decay (Bhattacharjee 2004) 


The primary questions addressed in this paper for incompressible MHD turbulence are also 
fundamentally important for the solar corona. In fact, our approach has many similarities with 


observational studies of solar flares (Crosby et al.||1993; Shimizu 

1995 

Boffetta et al.||1999 

Parnell 

Hi Jupp||2000 Hannah et al.||2008 

Aschwanden et al.||2000 ;|Uritsky et al.||2013, 2007 Veronig et al. 

2002; Aschwanden et al. 2014) and stellar flares (Benz Hi Giidel 

1994 

Audard et al. 1999 

Collura 

et al.|1988; Pallavicini et al.|1990 

Giidel et al.j2003 Telleschi et al.|2005 

). In these studies, the time- 


series of extreme UV, soft X-ray, and hard X-ray emissions from the Sun are used to characterize 
solar flares. Measured quantities include the size, duration, peak intensity, and fluence of the 
flares, from which the dissipated energy is inferred. In particular, the probability distribution for 
dissipated energy is of central importance due to its role in assessing the feasibility of the nanoflare 


model for coronal heating (Parker 1983 1988). This distribution is observed to obey a power law 


across eight orders of magnitude, with an index generally close to —1.8 (Aschwanden et al. 2000), 


although the precise value of the index varies significantly between different studies depending on 
the time period, region, type of emission, and methods used to identify the flares. Most importantly, 
this index is shallower than the critical value of —2, suggesting that nanoflares do not dominate 
the overall heating of the solar corona (Hudson 1991) ^ 


This analogy with the analysis of solar flares suggests a practical application for the statistical 
analysis of spatiotemporal structures: to assess whether turbulence plays a role in the intermittent 
energy dissipation of the solar corona. It is presently unknown whether the complex dynamics of the 
solar corona arise from self-organized criticality, turbulence, or some other phenomena ( jGeorgoulis 

. This problem can eventually be addressed by a better comparison 
of the various models with observations. Our methodology is ideal 
for making this comparison. As a hrst step in this direction, we compare the results in this paper to 


2005| lUritsky et al.||2007|[m3 


between numerical simulations 


^ Given the distribution for energy dissipation, P{E) ~ the critical index is derived by noting that the total 

energy dissipation, Etot oc EP{E)dE, scales with the lower bound i?min if a > 2 and with the upper bound 

^min 

-E'max if Ck <C 2. 
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the observed statistical properties of solar flares, collected from a number of previous studies. We 
find several nontrivial similarities in both cases, despite the anticipated differences between driven 
incompressible MHD turbulence and the essentially force-free dynamics of the solar corona. This 
suggests that MHD turbulence may play a role in the energetics of the corona, a possibility which 
should be investigated more carefully in future studies. 


The outline of this paper is as follows. We hrst describe the basic concepts and subtleties asso¬ 
ciated with the temporal analysis of structures in Section We then outline the technical aspects 
of our procedure, based on extending our previous framework for the statistical analysis of spatial 
structures ( [Zhdankin et al.||2013 2014), including detailed discussion of the algorithms involved, in 
Section In Section we introduce the measurements that are made on spatiotemporal dissipa¬ 
tive structures. The results of our analysis on the combined spatial and temporal intermittency of 
energy dissipation in 3D MHD turbulence are described in Section We discuss the implications 
of our results in Section which includes a comparison of our results with the observed statistical 
properties of solar flares. We summarize our conclusions in Section 


2. Background 

2.1. Intermittency in MHD turbulence 


MHD describes the macroscopic dynamics of an electrically conducting fluid, such as a plasma 
(see, e.g. Biskamp (2003)). Given a uniform background magnetic field Bq = Bqz that is strong 
relative to turbulent fluctuations (estimated by the rms value, b^ms), and assuming that gradients 
along the background held are small relative to those in the perpendicular direction, the incom¬ 
pressible MHD equations can be written in a reduced form as 


(I T • V||) . Vx) 

Vu 


-V^P + + ff 

0 , 


( 1 ) 


where = v ± b are the Elsasser variables (with directions strictly perpendicular to Bq), v is 
the huctuating plasma velocity, b is the fluctuating magnetic held (in units of the Alfven velocity, 
Va = Bq/ y/A-KpQ, where po is plasma density), P is the total pressure, V_l is the gradient in 
the (x,y) directions while Vy is the gradient in the ^ direction, and is the large-scale external 
forcing. For simplicity, the huid viscosity v is taken to be equal to the magnetic diffusivity rj, 
and both are assumed to be uniform in space. In the reduced MHD approximation, only the 2 - 
component of current density j = jz = z- \/±xb and vorticity lo = ujz = z ■ V ± x v are retained; 
these two scalar helds contain complete information to describe the dynamics as well as energetics. 

Energy loss from the system is governed by resistive dissipation and viscous dissipation, with 
respective energy dissipation rates per unit volume given by erj = rjp and Since the 

current density and vorticity are in direct correspondence with energy dissipation rates, they are 
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logical dynamical fields in which to study the intermittency of dissipation. In this work, we focus 
on resistive dissipation only for simplicity; it is straightforward to apply the methods to obtain 
similar statistical results for the viscous dissipation. 


Intermittency in strong MHD turbulence is characterized by the appearance of thin, quasi- 
2D ribbon-like structures in the current density, known as (electric) current sheets, with similar 
structures in the vorticity field. These current sheets are anisotropic in three directions, with the 
largest scale coinciding with the direction of the local mean field, which is the direction of the 
guide field to good approximation. Their lengths and widths are distributed mainly in the inertial 


range, while the thicknesses are strongly localized inside the dissipation range (Zhdankin et al. 


2013, 2014). The current sheets are notable for their role in heating, particle acceleration (Kowal 


et al. 2012), and magnetic reconnection (Biskam^|1986 ). 


2.2. General remarks for temporal analysis 

In this paper, we aim to study the temporal properties of intense, spatially-coherent structures 
produced by intermittency in 3D turbulence. For concreteness, we consider current sheets in MHD 
turbulence. Despite their prominence, intermittent current sheets occupy a small volume that 
represents the tail of the probability distribution for current density j{x). To systematically study 
these structures, we employ a framework based on identifying and characterizing the regions in 
space with current densities exceeding some fixed threshold value, jthr; which is the only parameter 
inherent to the methodology. Each structure is then represented as a 3D volume bounded by an 
isosurface of j. 

It is straightforward in principle to extend this procedure into the temporal realm by applying 
the same threshold criterion to the 4D spatiotemporal field j{x, t), thereby obtaining 4D spatiotem- 
poral structures. However, although simple in principle, it is challenging in practice to analyze a 
high-resolution 4D data set in this way. Assuming that the analysis is applied to post-processed 
data from well-resolved simulations stored on a computer system, the data storage conditions limit 
the total number of available snapshots, negatively affecting the possible time cadence of snapshots 
and size of the overall time interval. Hence, the time resolution of post-processed data must gen¬ 
erally be worse than the internal time resolution of the simulation. This is a serious issue since 
the data must be well-resolved in all dimensions in order to properly resolve and track structures 
across a wide range of scales. 

A related problem is that the primary memory limits the amount of data that can be loaded 
by the analysis program at any given time, requiring the numerical procedures to work with small 
pieces of the overall data set at a time. Therefore, a feasible temporal analysis should be based 
on first identifying the structures in the spatial dimensions of a given snapshot, and then tracking 
them through time, from their formation to their destruction. The algorithms required to perform 
this task are rather complex, and must be designed in a robust and efficient manner. In particular. 

























there needs to be an algorithm that accurately associates structures in one snapshot with their 
time-evolved counterparts in subsequent snapshots. A fundamental challenge here lies in the fact 
that there may not be a unique correspondence between a structure in one snapshot and a structure 
in the adjacent snapshot, due to mergers and divisions. 

The interactions between structures cause the notion of a coherent time-evolving structure to 
be ambiguous. Instead, the objects of central importance are the processes involving structures. 
These processes can be variously thought of as sets of interacting structures, as dissipative events, 
or as flares (if we equate the dissipated energy with outgoing radiation). More generally, they can 
be thought of as branched spatiotemporal structures. 


2.3. Classification of processes 


In order to build an intuition and facilitate the subsequent discussion of processes, we now 
describe a convenient classihcation scheme for processes. This allows processes to be visualized 
diagrammatically, which has some superficial similarities to Feynman diagrams for quantum me¬ 
chanical scattering processes ( Feynman|1948 ); however, none of the mathematical symmetries char¬ 
acterized by the Feynman rules carry over, since, to the best of our knowledge, isosurfaces in the 
current density described by the MHD equations contain no conserved quantities. Regardless, the 
following classification scheme is a simple way to describe processes and their complexity. 


We first introduce some terminology. We define a state to be an individual spatial structure 
at fixed time, which represents the basic building block of processes. We assume that the states 
are given at times spaced by an inhnitesimal increment dt. We also assume that there exists a 
map between all states at any time t to other states at time t — dt and t -|- dt, which represents the 
instantaneous temporal evolution of structures from one state to another state. We define a path 
segment to be a bijective (i.e., one-to-one) sequence of states under this map, which represents the 
coherent temporal evolution of an individual structure while it does not interact with any other 
structures. We then dehne a path to be a path segment with bijectivity breaking down only at 
the initial and final states in the sequence. A process is then described as a set of paths that are 
(non-bijectively) connected at their endpoints, which represents a set of interacting structures. 

Diagrams of some simple conceivable processes are shown in Fig. Here, we represent paths 
schematically as lines with an arrow marking the direction of time. If the path begins by sponta¬ 
neous formation, i.e., from a peak that grows to exceed the detection threshold, then we mark the 
beginning of the path with an O. If the path ends by spontaneous destruction, i.e., from a structure 
that recedes below the detection threshold, then we mark the end of the path with an X. The third 
possibility is for the path to start or end with an interaction. Interactions between structures are 
represented by vertices connecting sets of three (or more) paths. A process schematically consists 
of a set of paths connected by a set of vertices. 


The simplest process is the evolution of an isolated structure, i.e., a structure that is formed 



Isolated structure 
O-►- 




Division 



Merger 



Scattering 



Fig. 1.— Diagrams of some simple processes, where formation is represented by O, interaction by 
a vertex, and destruction by X. An isolated structure is a process with no vertices. Division and 
merger processes are the next simplest case, with a single vertex each. Higher-order processes such 
as loops and scatterings have a larger number of vertices or vertices with more paths. 
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and then destroyed without interactions. Isolated structures are described by a single path and 
have well-defined histories with well-defined properties. Other structures undergo at least one 
interaction. The two simplest processes involving an interaction are the division of one structure 
into two structures and the merger of two structures into one structure. Since an interaction is 
non-bijective, the structures in these processes do not have completely well-defined histories, so 
many of the quantities used to describe isolated structures are ambiguous. However, we will see 
that a meaningful set of more general characteristics can be introduced. 

We claim that the most logical approach for a temporal analysis is to study processes rather 
than individual spatially-coherent structures, i.e., paths, which lose their identity upon interacting. 
This is also the most conservative approach, as it requires no fundamental changes to the method¬ 
ology used for the statistical analysis of spatial structures at fixed time, and requires no ad-hoc 
assumptions to treat the interactions. We also find that the statistical trends are more robust for 
processes rather than paths. 


3. Methods 

3.1. Outline of procedure 

In this section, we describe our algorithms for the temporal analysis. Since the procedure in 
its entirety is rather complicated, we include a brief outline in this subsection. A more detailed 
description is presented in the rest of this section. The procedure rests on a hierarchy of steps: 
hrst, we find sets of contiguous points (above the threshold) in each snapshot to obtain the spatial 
structures; next, we hnd sets of bijectively-connected states to obtain the paths; finally, we find 
sets of connected paths to obtain the processes. A schematic of the final result is shown in Fig. 
where processes (colored in green) of varying complexity are identihed on the space-time lattice. 

1. Identify all states (i.e., spatial structures) in each snapshot, and represent tem¬ 
poral connectivity by constructing a map between states in adjacent snapshots. 

(a) Load initial snapshot and determine states (by using threshold algorithm). 

• Store constituent points of each state in a temporary array; 

• Perform measurements on states and store in a permanent array; 

(b) For k = 2,..., Nsna.p, do the following: 

• Load fth snapshot, determine states, and store constituent points and measurements; 

• Construct the evolution map, which associates states in snapshot k—1 with temporally- 
connected states in snapshot fc, and vice-versa (do this by comparing constituent 
points of structures in both snapshots); remove temporary array for snapshot A: — I; 


2. Identify paths. 
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(a) Obtain paths from states that form a bijective sequence under the evolution maps. 

(b) Perform measurements on paths by referencing the constituent states. 

(c) For each path, construct (non-bijective) map that identifies other paths connected to it. 

3. Identify processes. 

(a) Obtain processes from sets of connected paths. 

(b) Perform measurements on processes by referencing the constituent paths. 

(c) Treat processes that exist during the initial snapshot or final snapshot as special cases 
(incomplete processes). 


3.2. Identification of states (spatial strnctures) 


We first describe our algorithm to identify the states on each snapshot, i.e., the spatial struc¬ 
tures at fixed times. We define states as spatially-connected sets (i.e., clusters) of points with current 
density magnitudes larger than a fixed threshold, jthv- Two points on the lattice are considered 
spatially connected if the distance between the points is strictly less than two lattice spacings, i.e., 
one is contained in the other’s 26 nearest neighbors. 


We note that this is not the only conceivable way to define a state; one can, for example, 
use a variable threshold based on local field quantities such as the local peak current density 
( ]Zhdankin et ah' 2013). However, a fixed threshold appears to be the simplest approach, having 
only one free parameter (jthr) and no need for ad-hoc treatment of special cases (e.g., mergers). The 
main drawback of using the fixed threshold is that it generally gives a large number of unresolved 
structures due to peaks near the threshold; these do not exhibit physically meaningful scaling 
relations and must be carefully ignored. Fortunately, they are only manifest as noise in the low- 
intensity, small-scale regime of parameter space and generally have a negligible contribution to the 
total energy dissipation and volume of the structures. To some extent, unresolved structures are 
inevitable regardless of how they are defined, since there will always be a population of small, short¬ 
lived processes representing structures that barely cross the detection criteria. Filtering procedures 
can be applied to remove the population of unresolved, noisy structures residing near the threshold, 
but for simplicity we do not apply any filtering in our analysis. 


The definition of a state can affect how a given process appears; it may determine, for example, 
whether the close approach of two structures is registered as a merger or not. The threshold 
value may likewise have a similar effect. However, we would like to emphasize that the statistical 
conclusions should be, and in our experience are found to be, broadly consistent regardless of the 
method or threshold. 


Our algorithm identifies states by scanning the lattice for points with current densities above 
the threshold. For each such point found, the neighboring points satisfying |j| > jthr are identified, 
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then neighbors of those neighboring points satisfying \j\ > jthr, and so on, until no more points 
remain. The coordinates of the constituent points of the state are stored in an array for later use, 
and then the lattice is scanned for any additional states (while ignoring the points already found 
to belong to a state). Once every state in a snapshot is identified, measurements are performed on 
the states (see Section]^ and stored for later use. 


3.3. Temporal association between states in adjacent snapshots 

In this subsection, we consider the connectivity of states in the temporal dimension. We 
describe our algorithm for identifying the time-evolved counterparts of a state, i.e., the states 
in the adjacent snapshots that represent the evolved structure. The result is a (non-bijective) 
map associating the set of states in one snapshot with the set of states in an adjacent snapshot, 
representing the instantaneous temporal evolution of structures. 

Consider a given state in the A:th snapshot. To hnd the future counterparts of the given 
state, we iterate through all states in the subsequent snapshot and determine which ones have any 
constituent points that are spatially-connected to the given state. In other words, we look for states 
with points that coincide with or neighbor any of the given state’s points (using the 26 nearest 
neighbor criterion). Any such states are identified as future counterparts of the given state. 

Let the states in each snapshot be denoted by an index i G {1,2,..., A^fc}, where is the 
number of states in the snapshot. We construct an array that stores the indices of the corresponding 
future states in the {k + l)th snapshot, which we call the forward evolution map (i) where 
i = 1,2,, Nk. Thus, if the ith. state in snapshot k is associated with the jth state in snapshot 
k + 1, then T^{i) = j. If the zth state has no future counterparts, then we assign T^{i) = 0, 
representing the null state, which corresponds to the destruction of the structure. If the ith state is 
associated to multiple future states (i.e., it divides), then we assign (i) multiple values containing 
all of these future states. 

We perform a similar procedure to hnd the past counterparts of the given state. In this case, 
we iterate through all states in the preceding snapshot to determine which ones have points that 
are spatially-connected to the given state. Any such states are identihed as past counterparts of 
the given state. Likewise, we construct a backward evolution map, Tff{i), which identihes state i 
in the kih. snapshot with its corresponding past counterparts in the {k — l)th snapshot. If the ith 
state has no past counterpart, then we assign (i) = 0, which corresponds to the formation of 
the structure. If the ith state is associated to multiple past states (i.e., it results from a merger), 
then we assign (i) multiple values containing all of these past states. 

The entire set of evolution maps where k = 1,..., Wnap is constructed by iterating through 
the snapshots and applying the above procedure immediately after identifying the states in pairs 
of consecutive snapshots. After is constructed for all states in a given snapshot k, the array 
containing the constituent points of the states is deleted to free up memory before loading the next 
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snapshot. 


3.4. Identification of paths (tracking algorithm) 


In this subsection, we describe our algorithm for identifying paths across multiple snapshots, 
which is an important intermediate step before processes themselves can be identified. A path is 
abstractly defined as a sequence of states that is bijectively-connected under the evolution maps 
(see Sec. 2.3). From the algorithms in the preceding subsections, we have a sequence of Asnap 
snapshots (denoted by index fc), each with a set of states, along with the evolution maps 
associating the states in each snapshot with their time-evolved counterparts in adjacent snapshots. 
To identify paths, we must hnd the bijective sequences in {T^}- 


Consider the ith state in snapshot k. Bijectivity of the forward evolution map is satisfied if 
there exists a unique state j in snapshot A: -|- 1 such that T^{i) = j and T'^^(j) = i. Any two 
states satisfying this condition form a path segment. If this condition fails, however, then state i is 
the endpoint of a path. This happens either when the structure is destroyed (T^ {i) = 0) or when 
it interacts (either (i) or '^k+l (j) is multi-valued). Likewise, if there exists a unique state I in 
snapshot A; — 1 such that T^{i) = I and = i, then these two states form a path segment; 

otherwise state i is the beginning of a path due to either formation (T^ (z) = 0) or an interaction 
(either T^{i) or is multi-valued). Using these conditions, we can track the state of any 

given structure along a path through a sequence of snapshots, until an interaction is encountered. 


We iterate through the states in all Wnap snapshots and use the above tracking procedure 
to identify the associated path and its constituent states, marking those states so that they are 
ignored in the remaining iterations. As a result, we obtain a set of iVpath paths. For each path, 
we construct an array which contains the indices of the constituent states, which is referenced 
to perform measurements on the path. We also construct an array containing the indices of the 
predecessors, which are the other paths connecting to it from the beginning of the path. The 
predecessors are determined by operating with T~ on the first state of the path to obtain all of 
the past states, and finding the paths that contain these past states. In a separate array, we store 
the indices of successors, which are the other paths connecting to the end of the path, obtained 
by operating with T"*" on the final state of the path. The predecessors and successors of paths 
characterize the vertices between paths. Note that if the number of predecessors is zero, then 
the path is formed spontaneously. If the number of successors is zero, then the path is destroyed 
spontaneously. 


3.5. Identification of processes 

We finally describe how to identify processes from the set of paths and their predecessors and 
successors. Recall that processes are described as sets of connected paths. Therefore, we hrst 



- 14 - 


i k 

X 



Fig. 2.— A schematic of structures evolving in time, highlighted in green (shown in one-dimensional 
space for clarity). Structures from the initial and final states are marked in red. 



Fig. 3.— Schematic of structure evolution, shown in 2D space for clarity. The procedure stores 
the constituent points of the present state, checks the future snapshot for any states with points 
that are spatially-connected to the present state’s points, and then identifies these states as future 
counterparts of the present state. In the left panel, a structure evolves without interacting. In the 
right panel, a structure divides, having multiple future counterparts. This procedure is reversed in 
time to determine past counterparts, with a merger occuring for multiple past counterparts. 
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iterate through the set of paths. For each path, iterate through the predecessors and successors of 
the path, and then through the predecessors and successors of those, and so on, until no new paths 
can be obtained. The set of paths acquired in this way constitute a single process, and their indices 
are stored in an array corresponding to that process. The paths that have already been identified 
as belonging to a process are ignored in the remaining iterations of paths. 

Some processes will contain states from the initial snapshot or in the final snapshot of the 
dataset, which we call initial processes or final processes, respectively. These processes must be 
treated as special cases, since our information about them is incomplete. The simplest treatments of 
these processes are either to ignore them or to treat them as normal processes undergoing formation 
or destruction in the initial or final snapshots. For most of our analysis, we will ignore initial and 
hnal processes. However, they are included in the probability distributions for better statistics at 
long durations. Due to the relatively long interval of time in our simulations, the initial and final 
processes are a very minor contribution to the statistics, unless low thresholds are used. 

This concludes our discussion of the algorithms used to identify processes. We now have a 
sample of A^proc processes, each including an array of constituent paths. The paths contain all of 
the information necessary to perform measurements on the processes. These measurements are 
described in the next section. 


4. Measurements 
4.1. Measurements for states 

In this subsection, we describe measurements made on states, i.e., spatial structures at fixed 
time, which will be used in the next subsections to construct similar, more general quantities for 
paths and processes. 

The volume V of the structure is immediately obtained by counting the number of constituent 
points of the state and multiplying by the lattice volume element. The Ohmic energy dissipation 
rate is given by 

£ = J dVrjf , (2) 

where integration is performed across all constituent points of the given state. 

To characterize the geometry of each state, we measure the linear scales in three orthogonal 
directions. For length L, we take the maximum distance between any two points of the structure. 
For width W, we consider the plane orthogonal to the length and coinciding with the point of peak 
current density. We then take the maximum distance between any two constituent points in this 
plane to be the width. The direction for thickness T is then fixed by the orthogonality condition. 
We take the thickness to be the distance across the structure in this direction through the point of 



- 16 ~ 


peak current density. Since typical thicknesses tend to be comparable to the lattice spacing, we use 
a linear interpolation scheme to obtain finer measurements. Note that this method can cause the 
thickness to be over-estimated in some cases with complex morphologies, for example, in structures 
with S-shaped cross-sections, but these appear to be a relatively small population manifest as 
spurious high-value measurements. These definitions automatically satisfy L >W >T. 


We measure these scales in units of the system size in the direction perpendicular to the guide 
field. See Zhdankin et al. (2014) for more detailed information on this method applied to high- 
resolution simulations of MHD turbulence, as well as a discussion on alternative methods based on 
the Minkowski functionals. 


States can also be characterized by some discrete properties. One example is the direction of 
the current flow, i.e., orientation. Another example is the presence/absence of topological features 
such as X-points or 0-points [see Zhdankin et al. (2013) and Servidio et al. (2010)]. A final example 
is the Euler characteristic, i.e., genus. These characteristics will be relegated to a future paper. 


4.2. Measurements for paths 

In this subsection, we describe measurements for paths, which are conceptually simpler than 
those for processes due to the bijectivity condition. These measurements will be used and general¬ 
ized in the next subsection to characterize processes. 

The evolution of a path can be described by the time-series £{t), V{t), L{t), W{t), and T{t) 
of instantaneous characteristics defined in the previous section for the constituent states. Consider 
a path given by a sequence of states at times tk, where k G {1,..., W} and W is the number of 
constituent states. Assuming a fixed cadence of snapshots, the states are separated by a fixed time 
interval At. The A:th state has characteristics denoted by 14, L^, 114, and T^. 

One of the most basic properties of a path is its duration (or lifetime) r, defined as 

T = tN, -ti = {Ns - l)At. (3) 

The energy dissipation rate of states generalizes to the dissipated energy, 

. Ns 

E= dt£{t) = '^£kAt, (4) 

fc=i 

where time integration is performed over the duration of the path. We also define the peak volume. 
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peak energy dissipation rate, maximum length, maximum width, and maximum thickness as 

max(y(t)) = max({Vfc}) 
max(£i(t)) = max({<Sfc}) 
max(L(t)) = max({Lfc}) 
max(kk(t)) = max({H4}) 

max(T(t)) = max({Tfc}) (5) 

Note that since these quantities are local rather than time-integrated, they may be sensitive to 
chaotic fluctuations. As an alternative, we can consider time-averaged quantities, which however 
are less easily generalized for processes. 


kmax — 

f — 

^max — 

Amax — 

kbinax — 

T — 

^ max — 


4.3. Measurements for processes 

We now describe how to generalize the quantities defined for paths in the previous section 
to processes, i.e., sets of interacting structures. We characterize each process by the number of 
constituent paths, Np. Processes with a single path, Np = 1, are isolated structures. Processes 
with three paths, Np = 3, are division or merger processes. Processes consisting of more than three 
paths, Np > 3, are higher-order processes, containing either more than one vertex or vertices joining 
more than three paths. Other related measures of the complexity of a process include the number 
of vertices N^, the number of constituent states Ng, the number of internal paths Wnt (i-C- paths 
that begin and end in vertices), and the number of incoming paths A^in or outgoing paths Aout- 

Consider a process with constituent paths enumerated by index n = 1,..., Np, each extending 
from initial times tn to final times Let A„, Pmax,n, ^max,n, Wmax,n, and Lmax,n be the 

characteristics of the nth path as defined in the previous subsection. We define the process duration 
by 

T = max({4}) - min({t„}), (6) 

which is simply the time interval from the formation of the first existing constituent path to the 
destruction of the last existing constituent path. The total dissipated energy A of a process is de¬ 
fined by integrating rjj'^ across the enclosed 4D spacetime region. This is found in a straightforward 
manner from 

Np 

E = ^En. (7) 

n=l 

Likewise, we can define the peak volume Imax and peak energy dissipation rate Tmax as the maxi¬ 
mum of the peaks corresponding to the constituent paths, 

Lmax — aiax({ Vjnax,n}) 

^^max = max({Tmax,n}) • (8) 
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Note that an alternative definition for 14iax (^max) can be based on the maximum of the volume 
(energy dissipation rate) summed for all states belonging to the process at any given time. These 
two definitions may differ for processes with a large number of paths, but will otherwise give similar 
results; we use the first definition only for simplicity. The most difficult quantities to generalize for 
a process are the characteristic spatial scales. There appears to be no universally satisfactory way 
to obtain average characteristic scales for a general configuration with many paths. If we apply 
an average across all paths (or all states) constituting the process, then the result may be skewed 
toward unphysical short-lived paths (or states). A simpler alternative is to take the maximum scale 
corresponding to any path in the process, 

Tmax = Ill&xdZ/maXjn}) 
kbmax — Cnax({ VTniax,n}) 

^max — CCLax({rjnax,n}) • (9) 

One alternative definition for the spatial scales, which gives results consistent with the definition 
that we apply, is to take the spatial scales from the largest state at the moment of peak energy 
dissipation rate. 


5. Results 


5.1. Simulation details 


Before discussing the results of our temporal analysis of MHD turbulence, we first describe 
the numerical simulations. We consider simulations of strong, incompressible MHD turbulence 
driven at large scales. These simulations solve the reduced MHD equations, Eq. using a fully 


dealiased 3D pseudo-spectral algorithm (see Perez et al. 2012 for specific details on simulations). 
The ratio of magnetic guide field to rms fluctuations is set to Ro/^rms ~ 5. The periodic box is 
elongated in z (the direction of the guide field) by a factor of L||/L_l = 6, where L±_ = 27r is the 
perpendicular size of the domain in simulation units. Turbulence is driven at the largest scales 
by colliding Alfven modes, generated from statistically independent random forces in Fourier 
space at low wave-numbers 2'it/L± < kx,y < 2{2tt/ L±), kz = 2'k/L\\. The Fourier coefficients of 
are Gaussian random numbers with amplitudes chosen so that ~ 'I’rms ~ 1- The forcing is 
solenoidal in the perpendicular plane and has no component along Bq. The random values of the 
different Fourier components of the forces are refreshed independently on average about 10 times 
per eddy turnover time. The Reynolds number is defined by Re = Urms(T±/2vr)/i^; we set u = r] so 
that the magnetic Reynolds number, Rm = Vrms{L±/2Tr)/r], is equal to Re. Durations (and other 
timescales) are measured in terms of large-scale eddy turnover times of the turbulence, given by 
Teddy = L ±/~ 1- The analysis is performed on time intervals of durations rtot, ah of which 
begin after the simulations reach statistical steady state. 


We consider four simulations with parameters shown in Table Case 1 is a lower-resolution 
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(256^) run used to establish convergence of the results with resolution. Three independent runs 
with resolution 512^ are chosen to study scalings with Re, although the relatively limited range 
{Re = 800 — 1800) inhibits precise measurements of the scalings. Of these runs, Case 3 with 
Re = 1250 is the most robust data set, having the highest cadence, being nominally well-resolved 
dynamically, and having the longest time interval (rtot = 15.6). 

We analyze snapshots dumped at a cadence {At)~^, with At being larger than the internal 
time step in the simulation. For reference, we now establish the naive estimate for the minimum 
cadence required to properly track structures between two adjacent snapshots. This is estimated 
by requiring that the distance advected by the flow during At is less than the typical current 
sheet thickness. Estimating the former as v^msAt and the latter as l^rms/ithr; we require At < 
&rms/(^^rmsithr) ~ l/jthr- For the four cases in Table[^ we have jrms £ {12.1,11.8,14.6,17.4}, which 
gives the condition At < {1/12.1,1/11.8,1/14.6,1/17.4}jrms/jthr- The cadences given in Table 
fall short of satisfying this condition for thresholds more than three or four times larger than the 
rms fluctuations. However, we note that this condition is somewhat alleviated in practice because 
associations are made for states that do not fully overlap (since displacements of one lattice spacing 
satisfy the connectivity condition) and the condition is applied to all of the points in the (generally 
large) structure. Most of the results in our analysis show robust convergence with cadence or are 
only weakly sensitive. 


5.2. Global results 

We hrst describe the global features of energy dissipation in the simulations. In the hrst panel 
of Fig. we show the total ohmic energy dissipation rate in the system, Stot{t), during the given 
time intervals for all of the simulations in Table If we associate the dissipated energy with prompt 
optically-thin emission, i.e., if we assume that all dissipated energy is converted immediately into 
radiation in an optically-thin environment, then this represents a light curve for the system. The 
mean of £tot{t) is very close to 1.0, in agreement with the energy input from the large-scale forcing. 
The rms fluctuation about this mean is approximately 0.15 for all cases. In the second panel of 
Fig. 0 we show the power spectra of Stot{t), which exhibits a power law with index close to —2.0 
for all cases. 


Table 1: List of numerical simulations 


Case 

Resolution 

Re 

At 


1 

256^ 

800 

1/64 

10.0 

2 

512^ 

800 

1/32 

12.2 

3 

512^ 

1250 

1/64 

15.6 

4 

512^ 

1800 

1/32 

12.2 
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Fig. 4.— Left panel: the total energy dissipation rate in the system for the analyzed interval of 
time in the four simulations. Right panel: the power spectrum of this time series, showing a power 
law index near -2. The colors correspond to cases 1 (magneta), 2 (red), 3 (blue), and 4 (green) 
from Table m 




Fig. 5.— The fractions of total energy dissipation (left) and volnme (right) in structures with 
\j\ > jthr- The curves correspond to cases 1 (magenta), 2 (red), 3 (bine), and 4 (green) from 
Table Case 2 is consistent with the result from fully-resolved simulations obtained in [Zhdankm 

eTalldMIl ). 
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We now consider the energetics of structures with varying thresholds (normalized by the rms 
current density for the entire time interval, noting that jrms = \J£tot /"fiKot cx; The fraction 

of total energy dissipation occurring in structures with \ j\ > jthr, given by djP{j)rjj‘^ / (f^ djP{j)r]j 

where P{j) is the probably distribution of j, is shown in Fig. We find that this function has 
a wide tail characteristic of intermittency, declining exponentially at large jthr- We also show the 
fraction of total volume occupied by the same structures, given by djP{j), in the second panel 
of Fig. The fraction of energy dissipation always greatly exceeds the occupied volume, with, for 
example, about 30% of the energy dissipation occurring in 1% of the volume at jthr/jrms ~ 3. Our 
previous study has shown that the fractions of energy dissipation and occupied volume both may 
be universal for sufficiently well-resolved, high-i2e simulations ( Zhdankin et al.||2014 ). Cases 2 and 
3 closely match these previous results. Since significant deviations occur at large jthr/jrms for the 
other cases, they may not be completely well-resolved. 

Before proceeding to the quantitative analysis of processes in the simulations, we show an 
example of a relatively simple process in Fig. This process has a duration r ~ 0.5 with 31 
distinct paths. Snapshots of a few representative states (in green) are shown on a subdomain of the 
simulation lattice (with size 0.10 x 0.14 x 0.90). These images do not account for the elongation of 
the lattice along the (vertical) guide field, which would emphasize the ribbon-like character of the 
structures. The broadest part of the structure is shown from the given perspective. In addition 


to these snapshots, we also show the schematic diagram of the process (as established in Sec. 2.3), 


with the snapshot times marked in red. It is evident that this process is highlighted by a major 
division which occurs after the structure is stretched. A disproportionally large number of paths 
in this case are produced at the final stages of the process, as it decays toward the threshold. 


5.3. Aggregate quantities 

We now consider some aggregate quantities for the sample of processes in our simulations. 
Note that occurrence rates are generally not a very robust statistic, since they can be strongly 
affected by structures near the threshold, which in turn are strongly affected by resolution and 
cadence. A more robust analysis would filter out the small-scale and under-resolved structures to 
circumvent these numerical issues. We use no filtering in the present analysis, both for simplicity 
and because we are only interested in broad trends. 

Some general results are shown in Table for all of the simulations with fixed threshold 
ithr/jrms ~ 6.8 and cadence At~^ = 32. Here, Aproc is the number of processes, Apath is the 
number of paths, Aint is the number of processes with interactions (i.e., processes consisting of 
more than one path), Aisoi is the number of isolated structures (i.e., processes consisting of one 
path), Ajner is the number of merger events (where a vertex with n ingoing paths counts as n — 1 
mergers), A^iv is the number of division events (where a vertex with n outgoing paths counts as n— 1 
divisions), Afoi-m is the number of formation events (i.e., the number of paths with no predecessors), 
Ades is the number of destruction events (i.e., the number of paths with no successors), A 3 _vert is 
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Fig. 6.— An example of a typical process with duration r ~ 0.5, shown in green on the simulation 
lattice. A schematic diagram of the process is also shown, with the snapshot times marked by a 
red line. 


Table 2: Aggregate quantities in four simulations (ithr/jrms ~ 6.8, At = 1/32) 


Quantity 

256^, Re = 800 

512^, Re = 800 

512^, Re = 1250 

5123, iie ^ ;l800 

state) 

194 

288 

657 

1328 

-^proc 

914 

1271 

4272 

11608 

proc / Apath 

0.218 

0.240 

0.278 

0.339 

Aint/Aisol 

0.278 

0.329 

0.187 

0.115 

Amer / ^div 

0.839 

0.776 

0.800 

0.823 

Aform / A(Jes 

0.871 

0.853 

0.886 

0.911 

AQ-vert / A3-vert 

0.341 

0.353 

0.485 

0.523 
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the number of three-point vertices, and A^n-vert is the number of n-point vertices with n > 3. All 
of these preceding quantities are normalized to the number per eddy turnover time. We also show 
the mean number of states per snapshot, (iVgtate)- Table mainly shows the ratios of these various 
quantites; see the much more detailed Table in Appendix |A.1| for actual occurrence rates. 

We now make several remarks about Table First, (Wtate) and Aproc in Case 1 are smaller 
than the corresponding values for the higher-resolution Case 2, with a discrepancy of about 30%. 
In fact, there is a similar discrepancy in the number of occurrences for all measured quantities 
(see Table in Appendix). This can be attributed to the fact that Case 1 is somewhat under¬ 
resolved and therefore misses some of the small-scale dynamics present in Case 2. Second, (Agtate) 
and Aproc both strongly increase with Re, obeying estimated scalings of (Agtate) ~ Re^'^ and 
Aproc ~ Re^ '^. Third, three-point vertices are more common than higher-order vertices, although 
the higher-order vertices still occur in significant numbers. Higher-order vertices appear to be an 
unavoidable consequence of time discretization. Indeed, the ratio An-vert/As.vert increases with Re, 
implying that interactions occur over smaller timescales for higher Re. Fourth, although there are 
fewer processes with interactions than processes with no interactions, they contain the majority of 
the paths (i.e., Aproc/Apath < 1/2). Fifth, divisions are somewhat more common than mergers, 
implying a time asymmetry in the interactions. Equivalently, there are more destructions than 
formations. This asymmetry is reasonable since the time-reversal symmetry of the ideal MHD 
equations is broken by resistive and viscous dissipation. The preference for divisions over mergers 
may be a manifestation of the direct cascade of energy from large scales to small scales. 

We next consider the effect of cadence, At~^, on the above quantities. We show the various 
quantities from Case 3 for 1/64 < At < 1/4 in Tablej^ Most quantities monotonically increase when 
the cadence is increased, or equivalently, At is decreased (see Tablej^and Tablej^in Appendix |A.l[ ). 
However, Aproc fii'st increases then decreases with cadence, showing a local maximum near At « 
1/16. This unintuitive result may be explained as follows. Although the number of paths always 
increases with cadence, the connectivity of paths changes. For low cadence, paths tend to be 
isolated structures. This may be a sign that the cadence is insufficient to properly track structures 
- in particular, in the limit of very large At, snapshots become completely uncorrelated, so most 
processes appear as single states, therefore Apmc ~ (Ttot/At) {Nsmte) ■ For intermediate cadence 
(i.e., once structures are properly tracked), the paths interact more often when cadence increases, 
eventually decreasing Apmc due to the combining of isolated structures. The transition between 
these two regimes can be expected when Aproc/Apath ~ 1/2, which is close to the local peak in Aproc 
near At ~ 1/16 for this case. For very high cadence, Aproc either saturates or increases once again. 
The latter scenario may occur if additional isolated structures appear at the shortest time-scales, 
or if the processes are fractal. This trend is observed for Case 1 (see Table in Appendix |A.1[ ). 

Finally, we consider the effect of the threshold, jthr; on the above quantities. We show the 
various quantities from Case 3 for 5.5 < jthr/irms < 9.6 in Table For these relatively high 
thresholds, the occurrence rates for all quantities increase as jthr decreases, due to the larger 
sample size of paths (see Tablein Appendix|A.l[). In particular, we estimate that (Agtate) ~ /thr ^ 
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and A^path ~ itlir although these must deviate from a power law as jthr —^ jrms- In contrast, 
the ratios of the occurrence rates change relatively little. The asymmetry of interactions decreases 
and the relative proportion of isolated structures increases as jthr decreases, likely due to a larger 
sample of small structures near the threshold. 


The onset of percolation occurs near jthr/jrms ~ 6.8 for structures through space (in the 
periodic domain along the guide field, when max{L} —)• L||/(27r)), and near jthr/jrms ~ 5.5 for 
processes through the given time interval (when max{T} —>• rtot)- We call this latter quantity the 
percolation threshold, which, in the limit rtot —?• oo, is a fundamental characteristic current density 
of the system. For jthr below the percolation threshold, the initial processes and final processes 
(defined in Sec. 3.5) contain a large fraction of the dissipated energy, as can be seen in Table 
from F'interior/I^'aii) which is the ratio of energy dissipated by interior processes (i.e., processes that 
contain no states from the initial or final snapshots), FlinterioD to energy dissipated by all processes 
(including initial and final processes), Flail- This ratio is large (i.e., Flmterior/Flaii > 0.7) until the 
percolation threshold is approached near jthr/jrms ~ 5.5, where the ratio quickly becomes small 
(Flinterior/Flail ~ 0.29). The percolation threshold sets a practical limit on the smallest threshold for 
a reliable temporal analysis, since percolation otherwise interferes with the statistics of structures 
at the largest scales. 


5.4. Probability distributions and scaling relations 

We now describe the probability distributions for the process characteristics, defined in Sec. |4.3[ 
For clarity, we focus on Cases 2-4, which have resolution 512^ and varying Re. We choose a relatively 
high threshold of Jthr/jrms ~ 6.8, which is well above the percolation threshold. We retain initial 
and final processes for better statistics. The distributions are converged with respect to cadence 
and resolution (based on a comparison between Case 1 and Case 2 in Table [^, although it is 
possible that they have not fully converged with Re. 

We begin with the distribution for dissipated energy, P{E), shown in the first panel of Fig. 


Table 3; Variation of select aggregate quantities with cadence (Case 3: 512^, Re = 1250, jthr/jrms ~ 
6 . 8 ) 


Quantity 

At = 1/64 

At = 1/32 

At = 1/16 

At = 1/8 

At = 1/4 

^proc 

3311 

4272 

4908 

3704 

2197 

E proc/ E path 

0.120 

0.278 

0.562 

0.780 

0.895 

Wnt/Wsol 

0.429 

0.176 

0.067 

0.030 

0.016 

A^mer/-^div 

0.845 

0.799 

0.797 

0.771 

0.778 

-^form/-^des 

0.831 

0.886 

0.953 

0.978 

0.992 

^ n-vert 3-vert 

0.367 

0.485 

0.607 

0.697 

0.772 
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Table 4: Variation of select aggregate quantities with changing threshold (Case 3: 512^, Re = 1250, 
At = 1/64) 


Quantity 

Jthr/jrms ~ 9.6 

jthr/jrms ~ 8.2 

Jthr/jrms ~ 6.8 

Jthr/jrms ~ 5.5 

state) 

190 

343 

657 

1287 

^proc 

1105 

1777 

3311 

6314 

■^proc /-^path 

0.139 

0.125 

0.120 

0.116 


0.491 

0.481 

0.429 

0.383 

Nmer / Ndiv 

0.825 

0.832 

0.845 

0.863 

-^form/ -Ades 

0.820 

0.820 

0.831 

0.869 

n-vert /3-vert 

0.313 

0.331 

0.367 

0.400 

^'interior/ .Fall 

0.88 

0.86 

0.70 

0.29 

max {L} 

2.6 

3.3 

6.0 

6.0 

max{VF} 

0.32 

0.33 

0.44 

0.55 

max {r} 

3.3 

3.9 

8.5 

13.8 



Fig. 7. — The distributions for dissipated energy E, peak energy dissipation rate Tmax, and energy 
dissipation rate (of states) £. These have power laws with index close to —1.75 for P{E) and —2.0 
for T’(Tmax) and P{£)- The curves correspond to Re = 800 (red), Re = 1250 (blue), and Re = 1800 
(green). 
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We find that P{E) has a power law tail with index near —1.75 ±0.1. The power law region extends 
across approximately three orders of magnitude in E, from E ~ 10“® up to about E ~ 10“^. For 
smaller E, the distribution is shallower and appears to be non-universal, likely due to a combination 
of dissipation-range effects and threshold effects. With increasing Re, the power law extends to 
smaller E, consistent with a longer inertial range. If we instead consider the distribution for 
dissipated energy in isolated structures or in paths alone, then there is no clear power law. 


The distribution for the peak energy dissipation rate, T’(Tmax)) is shown in the second panel 
of Fig. We find that T’(Tmax) has a power law with index close to — 2.0 ± 0.1 in the range 
^max ~ 10“^ to Tmax ~ 10“^. Incidentally, this index is also observed in the distribution for energy 
dissipation rates, P{£), obtained from the population of spatial structures, i.e., states, shown in 
the third panel of Fig. The power law index of —2.0 for P{S) is in agreement with our previous 
analysis using much higher Re (Zhdankin et al. 2014). It can be shown that, if one assumes that 
all processes are single paths that evolve with identical, rescaled functional forms, then the indices 
for P(Tmax) and P{S) will in fact be equal. The assumption of identical evolution for processes 
of all durations is consistent with our results in Subsection 5.5, See Appendix A.2 for an analytic 
derivation of this relation. 


The distribution for process duration r is shown in the first panel of Fig. The durations 
extend to well above an eddy turnover time, sometimes comparable to rtot- The distribution from 
r ~ 0.2 to r ~ 8 can be fit by a power law with index near —3.2 ± 0.2. Likewise, the distributions 
for maximum length Lmax and maximum width VFmax have power laws with indices near —3.2, also 
shown in Fig. The distributions for all of these geometric quantities are related due to the strong 
correlations, described later in this subsection. 

We now remark on the distribution for maximum thickness. As shown in the first panel of 
Fig.§ F’f T'ma x) peaks at the lattice scale, h = 1/512 « 0.002 (in units of system size, L±). This is 
mainly due to the large population of under-resolved structures near the threshold. It is therefore 
more transparent to consider the distribution weighted by energy dissipation, i.e., to consider the 




Fig. 8.— The distributions for duration r, maximum length Lmax, and maximum width ITmax, all 
showing power laws with index near —3.2. The curves correspond to Re = 800 (red). Re = 1250 
(blue), and Re = 1800 (green). 
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total dissipated energy of processes at given Tmax, which we denote Etot{T^s.^). When -Etot(7"max) 
is compensated by Tmax, the maximum corresponds to the thickness scale at which most of the 
energy dissipation occurs. As shown in Fig F'tot(Fmax)7max is strongly peaked at T^ax a few 
times larger than the grid scale by (peaking at 2 to 3 times h, depending on Re). 


The extremely small values and variation of thicknesses makes it challenging to accurately 
infer the corresponding scalings in the given simulations. Higher resolution simulations are needed 
in order to make proper measurements. However, it is reasonable to expect that the bulk features 
of the structures and processes are insensitive to the thickness. Indeed, the robust scalings of the 
other characteristics are strong evidence in favor of this. This is tied to the fact that the other 
characteristics are essentially inertial-range quantities, which should be insensitive to the dissipation 
dynamics (assuming locality of the energy cascade). Further evidence for this comes from higher 
resolution simulations with similar Re, which show that the thickness of states is concentrated at 
a similar value as found here (Zhdankin et al. 2014). 


Finally, we consider the distribution for the number of paths per process, A^path; which is shown 
in Fig. We find that H(Apath) shows a robust power law with index near —2.0 ± 0.2. The most 
complex processes have ~ 10^ paths. It is remarkable that T’(A^path) shows such a robust power 
law across nearly the entire range of values (roughly 1 < Apath < 4 x 10^), since, a priori, it is not 
clear that the number of paths is a robust physical quantity. 


All of the distributions described above are insensitive to the threshold for large enough thresh¬ 


olds. As an example, we show P{£) and P{t) for 4.1 < jthr/jrms < 9-6 in Fig. 11 The distributions 
are similar in all cases with thresholds above the percolation threshold jthr/irms ~ 5.5. Deviations 
in the tails of both distributions are discernable when jthr/jrms = 5.5 and are more evident when 
ithr/jrms = 4:.l, Well below the percolation threshold. The percolation of processes steepens the 
tails of the distributions, consistent with undercounting the large-scale processes. 

We now describe the scaling relations between the various process characteristics. We show 


scatter plots of the different quantities versus process duration r in Fig. 12 For clarity, these are 
only shown for Case 3, with similar scalings for all other cases. We find that Lmax ~ SHVax ~ 0.6t, 


Srr 




^ r^, and A ~ (3 x 10“^)r^. We find that Tmax exhibits no evident correlation 
with r or other quantities. Therefore, to a good approximation, the thickness is constant. These 
scalings are then consistent with the simple geometric estimates, 


AjiiaxIFmaxFniax ~ T ; 

■f^max ~ I4iax??Ahr ~ ^ = I J dVr]f ~ based on the other correlations, 

assuming that the thickness and typical current densities are constant. 

A constraint between the indices of the distributions can be derived analytically if all pro¬ 
cesses are assumed to be single paths evolving with identical, rescaled functional forms. In this 
case, E ~ SmaxT is an exact relation; see Appendix A.2 for a derivation of this result. In addition, 
the distributions and scaling relations can be checked for self-consistency by using the conserva¬ 
tion of probability. For example, one may suppose that H(Tmax) ~ ‘^max) which implies that the 
distribution for peak energy dissipation rates is not dominated by weak or strong events. Then the 
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Fig. 9.— Left panel: the distribution for maximum thickness Tmax, which peaks at the lattice 
scale, h k. 0.02. Right panel: the same distribution weighted by dissipated energy and compensated 
by Tma x. showing that energy dissipation is dominated by processes with thicknesses a few times 
larger than the lattice scale. The curves correspond to Re = 800 (red), Re = 1250 (blue), and 
Re = 1800 (green). 



Fig. 10.— The probability distribution for the number of paths per process, which is well ht by a 
power law with index —2.0. The most complex processes contain thousands of paths. The curves 
correspond to Re = 800 (red). Re = 1250 (blue), and Re = 1800 (green). 
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measured scaling relations, £^max ~ and E ^ ^ iS’max, fix the indices of the other distributions. 

In this C8iS6, P(yE^ — -f^(‘^max) E which is relatively close to our measured index 

near —1.75. The scalings then also imply that P{t) ~ consistent with our measured index 
of —3.2. In general, it is clear that P{E) should be somewhat shallower than P{£max) due to 
integration across the duration, which increases with Tmax- 


5.5. Process evolution 


We now present results on the temporal evolution of individual processes. The following 
information is based on the time-series of instantaneous characteristics for each process, obtained 
from the constituent states at each snapshot. In general, the evolution of a given process is irregular 
and chaotic - in particular, long-lived processes are marked by frequent interactions and various 
phases of growth and decline. For example, the evolution of several characteristics for the two 
longest processes in Case 3 (at jthr/jrms = 6.8), which have durations of r 7.1 and 5.7, are shown 
in Fig. 13, These processes begin by rapid growth, followed by a relatively steady phase that is 


randomly kicked via interactions, and end by rapid decay toward the threshold. We investigate the 
evolution of a typical process by averaging over all processes of a given duration. 

To be concrete, we focus on the evolution of energy dissipation rate, £{t) for 0 < t < r. 


Shown in the left panel of Fig. 14 is the averaged energy dissipation rate normalized to peak energy 
dissipation rate, T(t/r)/Tmax) versus time normalized to duration, t/r, for processes of durations 
r G {0.125,0.25,0.5,0.75} in Case 3. The evolution is well approximated by a single sine mode, 
T(t/r)/Tmax ~ sin (vrt/r), independent of r. Since the long-lived processes have a similar evolution 
as short-lived processes (with time normalized to duration and energy dissipation rate normalized 
to the corresponding peak value), it is reasonable to average the statistics over structures with 
varying r. This type of average is shown in the right panel of Fig. 14 for all processes with r in the 
intervals {(0.25,0.5), (0.5,1), (1,2)}. The averaged T(t/r)/ Tma x continues to follow the same form 
up to r ~ 2, above which the statistical sample becomes limited. 

We next perform an average over processes of all durations to obtain (T(t/r)/Tmax), shown for 
Cases 2-4 in Fig. 15 It is clear that {£{t/ t) /Smax) ~ sin (vrt/r) holds to a very good approximation. 
For a more quantitative analysis, we show the power spectrum of {£{t/ t)/S^ asix) for Case 3 in the 
right panel of Fig. 15, which has a very steep decline in power going approximately as at low 


u, confirming that the w = 1 mode strongly dominates. The geometric characteristics V, L, W, 
and T show a similar temporal evolution as T, consistent with the strong correlations. 

Next we consider the temporal evolution of the instantaneous number of states involved in the 


process, Nstatesit/ t). This is shown for Case 3 in Fig. 16, with averages performed across durations 


in the four intervals {(0,0.5), (0.5,1), (1,1.5), (1.5, 2)} (left panel) and across all durations (right 
panel). The functional form of Ngtatesit) exhibits clear qualitative differences from £(t) and the 
geometric characteristics. The average across all durations can be approximated as a triangle 
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function, i.e., linearly increasing in time and then linearly decreasing in time. 


Although S{t) and Astates(i) are symmetric to a good approximation, a small asymmetric 
component can be discerned from Figs. 15 and 16 This asymmetry can be seen more clearly by 
decomposing the curve into symmetric and anti-symmetric parts. 


fsyrait) 


/asym(^) 


fji) + fjT - i) 

2 

fit) - fjr - t) 

2 


( 10 ) 


where f{t) is the given function on 0 < t < r. Taking f{t) = (T(t/T)/Tmax)) we show the 
symmetric and anti-symmetric parts in Fig. As noted before, the symmetric part is well fit by 
sin(7rt/r). We also find that it can be equally well fit by 1 — [{t — 0.5T)/(0.5r)]^'®, which is nearly 
indistinguishable from the sine peak. We hnd that the anti-symmetric part of (T(t/T)/Tmax) can 
be very well fit by sin(27rt/r), with an amplitude of 0.036. To investigate the asymmetry more 
precisely, we consider the first moments of the evolution curves, 


{t/'^)f 


Ioit/T)fit)dt 

fo fit)dt 


( 11 ) 


where deviation from 0.5 is indicative of temporal asymmetry. We show {t/T)£ and {t/ t) Ini' 
r < 1 in Fig. 18 We find that {t/T)£ is very close to but slightly below 0.5, while {t/r)Nutates is 
very close to but slightly above 0.5. At small r, the displacement from 0.5 initially grows with 
increasing r, but then asymptotes and becomes dominated by scatter at large r. Upon averaging 
over all durations, we find that for Re = {800,1250,1800}, {t/T)£ = {0.483,0.483,0.476} while 
{t/T)Nst^tes — {0.517,0.517,0.522}. Incidentally, the degree of asymmetry is comparable for both 
types of measurements, although in opposite directions. 


6. Discussion 


6.1. Comparison to solar flare observations 


In this subsection, we compare our statistical results for dissipative processes in MHD turbu¬ 
lence with the statistical properties of solar flares, taken from a number of observational studies. 
This comparison is not intended to draw any direct conclusions about solar flares, since our simu¬ 
lations are physically inadequate for describing the overall dynamics of the solar corona, although 
they may describe the turbulence that develops at small scales. Instead, the present comparison 
is motivated by the fact that solar flares are the best-observed natural example of intermittent 
energy dissipation in large-scale magnetized plasmas. For a more direct comparison, simulations of 
line-tied MHD ( Galsgaard fc Nordlund||1997 Ng &: Bhattacharjee 1998 Ng et al.|2012 Wan et al. 
2014) or other numerical models of the corona (Bingert Sz Peter 2011, [2^3) may be investigated. 
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The properties of solar flares are obtained from observations of extreme UV (EUV), soft X-ray, 
and hard X-ray emissions by applying a methodology similar to the one used in this paper. However, 
there are several important, unavoidable methodological differences that may affect the comparison. 
First, the emission may not be in direct association with the dissipation, making it nontrivial to 
infer the dissipation from the spectral amplitude. Indeed, although hard X-rays are thought to be 
promptly powered by dissipative magnetic reconnection events, soft X-rays and EUV can originate 
from aftereffects including chromospheric evaporation and cooling. Second, images of solar flares 
are projected onto a 2D plane, reducing the available information. In addition, there are several 
physical differences between driven, incompressible MHD turbulence and the solar corona. In 
contrast to volumetrically-driven turbulence in a periodic box, flares in the solar corona are generally 
modeled by force-free MHD with slowly-driven, line-tied boundary conditions. Additional plasma 
physics arising from kinetic and two-fluid effects may be important during magnetic reconnection. 
Following the reconnection event, other physical processes including chromospheric evaporation, 
radiative cooling, and thermal conduction may affect the decay of the solar flare. Nevertheless, we 
proceed with the comparison. 


In order to make a tangible comparison, we focus on a handful of studies which are method¬ 
ologically most similar to our present work. These are the papers by Uritsky et al. (2007|) (U07), 


Uritsky et al. (2013) (U13), and Aschwanden et al. (2014) (A14). These studies take extreme UV 


images of the corona and magnetograms of the photosphere to identify 3D spatiotemporal dissipa¬ 
tive processes. The range of indices for the distributions and scalings in these studies, along with 
our results, are shown in Table The statistics for dissipated energy and length agree favorably 
with our results, whereas they appear to differ for durations. The sub-diffusive growth of solar 


flares ( 

Aschwanden 

2012a; 

Aschwanden et al. 

2013 

), which has been modeled in the framework 

of self-organized criticality 

(Aschwanden 

2012b 

), also appears to be at odds with the evolution of 


processes measured in our work. 


Table 5: Comparison of distributions and scalings with solar flare studies 


Quantity 

MHD turbulence 

U07 

U13 

A14 

Index for P{E) 

1.75 

1.6- 1.7 

1.5 

1.8-2.2 

Index for P(Tmax) 

2.0 

- 

- 

2.1 - 2.5 

Index for P{L) 

3.2 

- 

2.5 - 2.9 

3.5-4.1 

Index for P{t) 

3.2 

1.9-2.1 

2.0 - 2.2 

2.2 - 2.6 

log E / log L 

3.0 

3.0-3.6 

3.0-3.1 

2.5 - 2.6 

log r/log L 

1.0 

1.8-2.3 

1.2 - 1.4 

- 


Comparing to more general studies of solar flares, our distribution for dissipated energy, with 
index near —1.75 ± 0.1, is close to the analogous measurements for total energy released in solar 
flares identified from hard X-rays, generally having an index quoted to be between —1.7 and —1.8 
(e.g., Aschwanden et al.|2000 Bromund et al.|[l995 Christe et al.|[2008). Similarly, our distribution 
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for peak energy dissipation rate fmax with index —2.0 ± 0.1 is close to observations of peak hard 


X-ray flux (e.g., Bromund et al. 1995) and soft X-ray flux (e.g., Aschwanden Sz Freeland 2012). 
Our distribution for duration with index near —3.2 ± 0.2 is somewhat steeper than the indices 
ranging between —2.2 and —3.0 for solar flare durations (Crosby et al.|[T9M Bromund et al.|[l995 


Veronig et al. 2002), although it is closer to the index for rise times, given as —3.4 in Christe 


et al. (2008) and —3.2 (during solar maximum; shallower otherwise) in Aschwanden & Freeland 


(2012). One final point of comparison is the asymmetry of processes, recalling that in our case, a 


process tends to dissipate slightly more energy at early times than late times. This is qualitatively 
in agreement with observations of solar and stellar flares, although the asymmetry appears to be 
much more pronounced in flares. The asymmetry can be dehned from the rise time trise and decay 
time tdecay as Aev = (tdecay “ Uise) / (tdecay + trise)] this is found to be 0.2 for X-ray flares, giving a 
peak at approximately 40% of the flare duration ( Christe et al.|[^08 ). In contrast, the asymmetry 
of processes in our simulations, based on this dehnition, is 0.034. 

In summary, the energetic and geometric statistical properties of dissipative events in MHD 
turbulence are consistent with solar flare observations, whereas the durations and temporal asym¬ 
metries present a noticeable discrepancy. The differences may be due to the Neupert effect, in 
which the chromospheric evaporation prolongs the decay of a flare observed in soft X-rays relative 


to hard X-rays (Neupert 1968; Dennis Sz Zarro 1993). This would explain why the distribution of 


process durations in our present work matches the distribution of solar flare rise times better than 
their total durations. 

The nontrivial similarities in the statistical properties between dissipative events in MHD tur¬ 
bulence and in the solar corona leaves open the possibility that MHD turbulence plays a governing 
role in the intermittency of the coronal energetics. This possibility has been advocated in numerous 


past studies (e.g., Georgoulis 2005 Uritsky et al. 2007, 2013) as an alternative to self-organized 


criticality. A more careful three-way comparison of the temporal statistics of dissipative events in 
MHD turbulence, self-organized criticality, and observations of the solar corona is left for future 
consideration. 


6.2. Implications for MHD turbulence 

The methodology developed in this paper is complementary to the tools conventionally used 
to probe turbulence and its intermittency. It also has several advantages over the other methods, 
despite the relatively complex numerical implementation and the necessity of a large data set. For 
example, given a relatively meager Re = 800, the distribution for dissipated energy per process 
shows a power law across nearly three decades in energy, whereas the inertial range is barely dis- 
cernable in the corresponding energy spectrum at the same Re. This large separation of scales may 
be attributed to the additional information given by the temporal dynamics of the turbulence. The 
analysis of structures also naturally describes the anisotropy and inhomogeneity of the dynamics, 
which is often challenging for other methods. 
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In this work and in our previous study of spatial structures (Zhdankin et al. 2014), we found 


that the distribution for the energy dissipation rates of states has an index close to the critical 
index of —2. This suggests that, at any given time, structures of all energy dissipation rates (in 
the inertial range) contribute almost equally to the overall energy dissipation rate. This may be 
a manifestation of the scale-invariance of inertial-range turbulence, since a distribution for energy 
dissipation rates with the critical index is equivalent to the energy dissipation being evenly spread 


across structures of all lengths (Zhdankin et al. 2014). 


Our present work goes beyond this previous result by establishing that the dissipated energy in 
evolving structures (processes) has a power law distribution with index shallower than the critical 
index, namely, with an index near —1.75 ± 0.1. This implies that intense dissipative events, i.e., 
large-scale and long-lived coherent structures, dominate the overall energy dissipation. This is a 
consequence of the linear scaling of duration with maximum length, which causes the distribution 
of dissipated energy to be shallower than the distribution of energy dissipation rates. 

The distributions, scalings, and evolution of processes appear to be insensitive to the Re sam¬ 
pled in our simulations. This suggests that these statistics may be valid for asymptotically large Re, 
relevant for space and astrophysical turbulence. This also implies that the results may be universal, 
i.e., insensitive to the mechanisms of energy input and dissipation, although this should be verified 
in the future by varying the boundary conditions, forcing mechanisms, and dissipation mechanisms. 
Examples of astrophysically-relevant driving mechanisms include the magnetorotational instability 
for accretion disks ( ]Balbus fc Hawl^ 1991| and line-tied driving for the force-free solar corona. 
It is also possible for the nature of intermittency to undergo a transition at sufficiently large Re, 
due to instabilities for large and morphologically complex structures. Therefore, it is important 
to verify our results in future simulations of MHD turbulence with larger Re, where more precise 
power-law fits can be obtained and a systematic study of the Re dependence can be investigated. 
It is challenging to apply our methodology, in its present form, to direct numerical simulations 
with larger Re, since both the spatial resolution and time cadence must be increased, making it 
impractical to store the full sequence of data snapshots. It may be necessary to perform the bulk of 
the analysis in parallel with the simulations, rather than analyzing post-processed snapshots as was 
done in this work. Alternatively, the amount of information used for the analysis may be reduced 
by, e.g., filtering out large-scale or small-scale modes. This is left for future consideration. 

Our methodology provides a new avenue to investigating temporal asymmetry, which was 
previously inferred in studies of MHD turbulence through the third-order moment or rate of energy 
flux (e.g., MacBride et al.||2008 Podesta||2008 Wan et ^|2010 ) and field-line diffusion (Beresnyak 


2014). The temporal asymmetry in our case is measured in the larger number of divisions than 


mergers, the tendency of the number of states in a process to be larger at late times than early times, 
and the tendency of the energy dissipation rate and geometric characteristics to be larger at early 
times than late times. Temporal asymmetry may occur from the onset of an instability, such as the 
tearing instability in resistive MHD or avalanches in critically self-organized systems. However, the 
asymmetry measured in our simulations is relatively small and does not significantly increase for 
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larger structures that may cross an instability threshold. Furthermore, the processes do not exhibit 
the impulsiveness expected when an instability is triggered. Therefore, we find it unlikely that the 
instability of structures plays a role in our simulations, although it will be an important signature 
to search for in future studies. Indeed, the tearing instability is expected to occur for laminar, 2D 
current sheets when they become sufficiently thin, which may occur at Rm ~ 10^ (e.g., Loureiro 


et al. 2007 Bhattacharjee et al. 2009 Uzdensky et al. [2010 ). It is conceivable that current sheets 


in a turbulent medium become unstable at different (possibly lower) Rm than naively expected 
([Loureiro et al.|2009); it is also possible that the instability is entirely absent. 


In our case, the temporal asymmetry may be linked to the turbulent energy cascade. Specif¬ 
ically, the inertial range of 3D MHD turbulence is characterized by a direct energy cascade from 
large scales to small scales. Therefore, turbulent eddies cause large structures with inertial-range 
lengths and widths to be broken into smaller structures, leading to a surplus of divisions over 
mergers, as well as more states at late times in a process. This can also explain why the energy 
dissipation rate and geometric characteristics are larger at early times, since a single large state 
may accomodate a higher current density than many individual states. Since the dynamics are oth¬ 
erwise time-symmetric in the inertial range, this asymmetry can be relatively weak. We note that 
another distinct contribution to the asymmetry can be from the dissipative term directly (rather 
than the cascade through the inertial range), relevant for states with lengths and widths that are 
near the dissipation scale. It is left to future work to better quantify the asymmetry and its origins, 
and to relate the measured quantities to the energy cascade rate. 


7. Conclusion 

The methodology developed in this paper, partly inspired by similar approaches in observa¬ 
tional astrophysics and self-organized criticality, presents a lucid picture of intermittency in MHD 
turbulence with dissipative structures playing a central role. Previously, this picture was incom¬ 
plete, based on a fraction of the available information. By exploring the temporal dimension, we 
completed this picture and found a richer, more valuable perspective of turbulence. 

This work demonstrates that the statistical analysis of spatiotemporal dissipative structures, 
i.e., processes or flares, can lead to concrete physical insights for intermittency in turbulence. The 
methodology is a natural but nontrivial extension of the methodology applied to characterize spatial 
dissipative structures in previous studies. We applied this analysis to numerical simulations of 
MHD turbulence to arrive at the following basic conclusions, which are insensitive to the numerical 
parameters (cadence, resolution, and threshold). We found that resistive energy dissipation occurs 
in current sheets that participate in intense, complex, long-lived processes with durations that 
may span several large eddy turnover times. These processes are analogous to flares in the solar 
corona and other astrophysical systems. The durations of the processes are directly proportional 
to the maximum lengths, providing a strong link between the spatial and temporal behavior. The 
energy dissipated in these intense processes is distributed as a power law with index near —1.75, 
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implying the dominance of large, intense flares. Incidentally, this index is consistent with observed 
energy distributions of solar flares. Processes are weakly asymmetric in time, dividing more often 
than merging. The averaged temporal evolution for the energy dissipation rate (and geometric 
properties) of the processes exhibits a nearly time-symmetric, sine-like form that is applicable to 
processes of all durations that were robustly sampled in our study. 

The insights gained from this methodology may help to understand magnetic reconnection, 
heating, and particle acceleration in space and astrophysical plasmas. Eventually, this methodol¬ 
ogy may also help narrow the imposing gap between our theoretical knowledge of energy dissipation 
in magnetized plasmas and the observational problem of coronal heating, and to uncover the sim¬ 
ilarities and differences between self-organized criticality and turbulence. In the future, we hope 
that the additional utilization of the methods described in this work will further evolve our present 
picture of turbulence. Indeed, the picture presented in this paper is incomplete in that it sidesteps 
the intermittency of vorticity and viscous dissipation, as well as its coupling to the current density 
and resistive dissipation, which should be scrutinized in a later work. In addition, we believe that 
further phenomenological modeling of the spatiotemporal dissipative structures in MHD turbulence 
is a promising pursuit. 
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A. Appendix 

A.l. Tables of occurrence rates 

The tables in this Appendix show the occurrence rates used to compute ratios in Tables EE 
and 1^ Table compares all cases, Tables and show varying cadence for Case 1 and Case 3, 
and Table shows varying threshold for Case 3. 


Table 6: Aggregate quantities in four simulations (ithr/jrms ~ 6.8, At = 1/32) 


Quantity (per Teddy) 

256^, Re = 800 

512^, Re = 800 

512^, Re = 1250 

512^, Re = 1800 

Apath 

4202 

5295 

15392 

34275 

^ proc 

914 

1271 

4272 

11608 

Aisol 

715 

956 

3600 

10416 

Aint 

199 

315 

672 

1193 

Adiv 

1494 

1663 

5352 

11585 

^ mer 

1253 

1290 

4276 

9540 

Ades 

1746 

2458 

7088 

17312 

Aform 

1520 

2096 

6283 

15766 

A3-vert 

1449 

1527 

3886 

7745 

^ n-vert 

494 

539 

1884 

4051 

(-^state) 

194 

288 

657 

1328 


Table 7: Aggregate quantities with changing cadence (Case 1: 256^, i2e = 800, jthr/jrms ~ 6.8) 


Quantity 

At = 1/64 

At = 1/32 

At = 1/16 

At = 1/8 

At = 1/4 

■^path 

6852 

4202 

2426 

1360 

720 

^proc 

1136 

914 

959 

921 

629 

Aisol 

894 

715 

834 

867 

613 

Aint 

241 

199 

126 

54 

16 

Adiv 

2614 

1494 

666 

199 

39 

^mer 

2277 

1253 

514 

141 

27 

Ades 

2295 

1746 

1418 

1078 

663 

Aform 

1984 

1520 

1276 

1022 

651 

As-vert 

2953 

1449 

512 

125 

24 

^ n-vert 

775 

494 

235 

71 

14 

state) 

194 

194 

193 

192 

189 
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Fig. 11.— The probability distributions P{S) (left) and P{t) (right) at various jthr/jrms- The 
distributions are affected by percolation through the time interval for thresholds below jthr/jrms = 
5.5 (green), visible in the curve for jthr/irms = 4.1 (magenta). 


Table 8; Aggregate quantities with changing cadence (Case 3: 512^, Re = 1250, jthr/jrms ~ 6.8) 


Quantity 

At = 1/64 

At = 1/32 

At = 1/16 

At = 1/8 

At = 1/4 

■^path 

27500 

15390 

8731 

4748 

2454 

^proc 

3311 

4272 

4908 

3704 

2197 

A'isol 

2316 

3600 

4600 

3595 

2163 

iVint 

994 

672 

308 

109 

34 

A'div 

11900 

5352 

1816 

495 

116 

^tner 

10060 

4276 

1447 

382 

90 

Ades 

8364 

7088 

5964 

4031 

2283 

Aform 

6954 

6283 

5685 

3942 

2264 

As. vert 

10750 

3886 

1099 

261 

57 

^ Q-vert 

3949 

1884 

667 

182 

44 
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Fig. 12.— Scatter plots of maximum length Fmax, maximum width VFmax, peak energy dissipation 
rate fmax, peak volume 14iax (relative to the system volume), and dissipated energy E versus the 
process duration r. 
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Fig. 13.— The evolution of several characteristics for the two longest processes for the Re = 1250 
case. The curves correspond to energy dissipation rate £ (black), volume V (blue), length L (red), 
width W (green), and thickness T (magenta). 




Fig. 14.— Energy dissipation rate versus time, averaged for all processes of given durations (left) 
and all processes in given intervals of durations (right). Also shown in black is the fit by a sine 
function. 































- 40 - 




Fig. 15.— Left panel: the evolution of energy dissipation rate {£{t/ t)/S max) versus time t/r, with 
the average performed across processes of all durations. The fit by sin(7rt/T) (shown in black) 
works very well. The colors correspond to Re = 800 (red), Re = 1250 (blue), and Re = 1800 
(green). Right panel: the power spectrum of {£{t/ t)/S^ nax) for the Re = 1250 case, showing a very 
steep descent as a power law with index near —5.0 at low ui. 




Fig. 16.— The instantaneous number of states in the process, Nstates{t/ t), across the duration of 
the process, averaged for processes with durations in given intervals (left) and for processes with 
all durations (right). 
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Fig. 17.— Left panel: the symmetric part of (<?(t/T)/<5max)) averaged across processes with r < 1. 
The fit by sin (irt/T) (shown in black) and 1 — [(f — 0.5r)/(0.5r)]^'® (shown in red) both work very 
well. Right panel: the corresponding antisymmetric part, fit by 0.036 sin (27rt/r) (shown in black). 




Fig. 18.— The first moment, (f/r), of the evolution of energy dissipation rate £{t/T)/ (left) 
and number of states Nsta,tes{t/ t) (right), versus process duration r for Re = 800 (red). Re = 1250 
(blue), and Re = 1800 (green). The curves are smoothed for clarity. 
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A.2. Constraints between indices of distributions 


If we suppose that all processes are single paths that evolve with identical (but rescaled) func¬ 
tional forms for the energy dissipation rate, then we can derive nontrivial constraints between the 
indices of the various distributions given in Sec. 5A Specifically, for simplicity, assume that all pro¬ 
cesses of duration r consist of a single path with the energy dissipation rate Sr{t) = £ma.x{T) f {t/ t) 
for 0 < t < r, where the universal shape function f{x) satisfies /(O) = /(I) = 0, 0 < f{x) < 1 and 
sup{/(x)} = 1 for 0 < X < 1. If we assume power-law distributions for all quantities, then they 
are given by 


P{E) ~ E-°‘ 

P{S) ~ £-^ 

P{£m!^x) ~ 

P{t) ~ T~'^ , (AI) 

with the constraint 

(/3-l)(7-a) = ( 7 -l)(a-l). (A2) 

The corresponding scaling relations are given by 

fma. ~ . (A3) 

Normalizable distributions require 7 > a. For the measured value of /3 = 2, 7 = 1/(2 — a) is 
required, so that 1 < a < 2. An example set of parameters satisfying this constraint are /3 = 2, 


Table 9: Aggregate quantities with changing threshold (Case 3: 512^, i?e = 1250) 


Quantity 

/thr/jrms ~ 9.6 

jthx/ Jrms ~ 8-2 

/thr/jrms ~ 6.8 

/thr/jrms ~ 5.5 

-^path 

7975 

14270 

27500 

54630 

^ proc 

1105 

1777 

3311 

6314 

A'isol 

694 

1200 

2316 

4565 

iVint 

341 

577 

994 

1749 

A^div 

3077 

5757 

11900 

26390 

^ mer 

2539 

4788 

10060 

22785 

Ndes 

2694 

4705 

8364 

12540 

A^form 

2210 

3858 

6954 

10900 

A^ 3-vert 

3081 

5557 

10750 

22120 

^ n-vert 

964 

1837 

3949 

8854 

state) 

190 

343 

657 

1287 
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7 = 3, and a = 5/3, along with scalings E ^ and £^max ~ all of which are consistent with 
the results in this paper. 

The derivation is as follows. We first relate the distribution of instantaneous energy dissipation 
rates, P{£), measured from states at random, to the distribution of peak energy dissipation rates, 
-P(Tmax); measured from processes. Assuming that one samples a random value £ from £r{t) = 
£max.{T)f {t/ t) with Uniform time sampling, so P{t) = 1/r, we obtain the distribution of energy 
dissipation rates from a process of duration r. 


P{£\t) 


dt 


d£ 


Pit) 


1 




E 


dx 

dfix) 


1 f <f’max('r) \ 

<?max(T)^V ^ ) 


(A4) 


where Xi [i = 1,... ,n) are the n roots of /(xj) — <?/<?max(x), and we have defined the function giy). 
The total distribution of energy dissipation rates is then 


Pi£) = 


/ OO 

dTPiT)P{£\T) 

min 


' min 
roo 


(irP(r) 
d£ 


£m 

Pi£n 


‘f^max(x ) 

£ 


£n 


On 


£ 


(A5) 


where Tm\n is defined such that £^max(xmin) = £■ The lower bound of the integral is required since 
processes with durations r < Tmin do not reach high enough ener gy d issipation rates to contribute 
to the distribution. Now assume that P(£^ma} 


£-d 

omax 


Then Eq. 


A5 


becomes 


Pi£) ~ d£m..£-i-^g 

/ OO 

dyy~^~^giy) 

~ £-d, (A6) 


where y = £raa.x/£- Therefore, assuming the integral in Eq. A6 converges, P(<?) has the same index 
as P(<f ma x). To relate this to other indices, note that the dissipated energy E(r) per process is 
given by 


dt£rit) = £^max('r)r [ dxf (x) ~ £^max('r)r , (A7) 

Jo 

where x = t/r, and the integral dxf (x) evaluates to a constant of order unity; hence, E ~ £m»^ T 
is exactly satisfied. Assuming P(r) ~ and P(E) ~ E~°‘, we can find the exponent A for E ~ 



— = ^ t-A-1 = ^ ^ ^ A = ^ 

dr P{E) T—\a a — 1 


(A8) 
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Hence, 8^^.^ ~ E/t ~ ^ ^( 7 -«)/((^-i) and 

C^<tmax 

Therefore we have /3 = 1 + (7 — l)(a — l )/(7 — a), as required. 
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